home *** CD-ROM | disk | FTP | other *** search
/ Power Bytes: Money & Finance / PowerBytes Money and Finance CD-ROM 01 / PowerBytes Money and Finance CD-ROM 01.iso / MacII / Accuracy / Accuracy.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1988-05-23  |  15.2 KB  |  566 lines  |  [TEXT/TPAS]
  1. { ACCURACY.PAS: 8087 math accuracy test, v 1.7, Oct 1987
  2.     (c) 1987, PC Tech Journal and Ziff Communications Co.
  3.                     written by Jim Roberts.
  4.  
  5. Uncomment the correct constant declarations for your compiler.
  6. Insert correct version number in `compil' variable,
  7. name of hardware system & math coprocessor in `machin' variable.
  8.  
  9. }
  10. { $MC68020+}
  11. { $MC68881+}
  12. {$R+}
  13. {$U-}
  14.  
  15. program ACCURACY2;
  16.  
  17.  
  18. USES 
  19.    
  20.       Memtypes, QuickDraw, OSIntf, ToolIntf, PackIntf, SANE, MacExtras;
  21.  
  22. TYPE
  23.    accarray  =  array[1..5] of extended ;
  24.    
  25.  
  26. CONST
  27.    compil   =  'Macintosh Turbo Pascal' ;
  28.    machin   =  'Macintosh II'  ;
  29.    minerr   =  1.0E-17 ;
  30.    logmin   =  17.0 ;
  31.    n        =  10 ;
  32.    y        =  1.0 ;
  33.    step     =  0.2 ;
  34.    iter     =  20 ;
  35.    itertrig =  5 ;
  36.    log10e   =  0.43429448190325182765 ;
  37.    pi       =  3.14159265358979323846 ;
  38.    piO2     =  1.57079632679489661923 ;
  39.    root2    =  1.4142135623730950488 ;
  40.    root3    =  1.7320508075688772935 ;
  41.    sqrtO2   =  0.7071067811865475244 ;
  42.  
  43. VAR
  44.    i, j, k, l, m, ntest,z      : integer ;
  45.    a,b,c                       : array[1..N,1..N] of extended;
  46.    sum, X                      : extended ;
  47.    xx, zz, quot                : extended ;
  48.    a0, a1, d0, d1, frac        : extended ;
  49.    p, p2                       : extended ;
  50.    th, err, logerr, diverr,
  51.      val, funct                : accarray ;
  52.    testerr                     : array[1..10] of extended ;
  53.    toterr                      : extended ;
  54.    done                        : Boolean;
  55.    s                           : DecStr;
  56.    textWidth                   : INTEGER;     { Maximum width of a character }
  57.    textHeight                  : INTEGER;     { Height of a character }
  58.    f                           : DecForm;
  59.  
  60.  
  61.  
  62. PROCEDURE DrawAt( x, y : INTEGER; s : Str255 );
  63.  
  64.    { Draw string at this coordinate, similar to an x,y location }
  65.    { on a conventional computer terminal. }
  66.  
  67.       CONST
  68.       
  69.          Xoffset  = 5;     { Blank pixels in left border }
  70.          Yoffset  = 12;    { Blank pixels in top border }
  71.  
  72.       BEGIN
  73.          MoveTo( Xoffset + x * textWidth , Yoffset + y * textHeight );
  74.          DrawString( s )
  75.       END; { DrawAt }
  76.       
  77. PROCEDURE NewFont( fontNumber, pointSize : INTEGER );
  78.  
  79.    { Select new font and adjust variables for DrawAt }
  80.  
  81.  
  82.       BEGIN
  83.          TextFont( fontNumber );          { Use new font }
  84.          TextSize( pointSize );           { Use new size }         
  85.          textHeight := pointSize + 2;     { Calc new character height }
  86.          textWidth := CharWidth( 'M' )    { Calc width of widest char }      
  87.       END; { NewFont }
  88.       
  89.  
  90.  
  91. procedure filla;
  92.    var
  93.       i,j     : integer;
  94. begin
  95.    for i := 1 to N do
  96.       for j := 1 to N do
  97.          if i <> j then a[i,j] := Y
  98.                    else a[i,j] := X + Y;
  99. end;
  100.  
  101. procedure fillb;
  102.    var
  103.       i,j     : integer;
  104.       f,d     : extended   ;
  105. begin
  106.    f := X + N*Y ;
  107.    d := 1.0/(X*f)     ;
  108.    for i := 1 to N do
  109.       for j := 1 to N do
  110.          if i <> j then b[i,j] := -Y*d
  111.                    else b[i,j] := (-Y+f)*d    ;
  112. end;
  113.  
  114. procedure fillc;
  115.    var
  116.       i,j     : integer;
  117. begin
  118.    for i := 1 to N do
  119.       for j := 1 to N do
  120.          c[i,j] := 0.0;
  121. end;
  122.  
  123. procedure mult;
  124.    var
  125.       i, j, k         : integer;
  126. begin
  127.    for i := 1 to N do
  128.       for j := 1 to N do
  129.          begin
  130.          sum := 0.0;
  131.          for k := 1 to N do
  132.             sum := sum + a[i,k] * b[k,j];
  133.          c[i,j] := sum;
  134.          end;
  135. end;
  136.  
  137. procedure sumit ;
  138.    var
  139.       i, j            : integer;
  140. begin
  141.    sum := 0.0;
  142.    for i := 1 to N do c[i,i] := c[i,i] - 1.0 ;
  143.    for i := 1 to N do
  144.       for j := 1 to N do
  145.          sum := sum + abs(c[i,j]);
  146. end;
  147.  
  148. function osgn(n:integer) : extended ; {negative if n odd}
  149. begin
  150.    if n = 2*(n div 2) then osgn := 1.0
  151.                       else osgn := -1.0 ;
  152. end;
  153.  
  154.  
  155.  
  156. procedure arith ;
  157.    {TEST 1: well-conditioned combinatorial matrix times its inverse.}
  158. var
  159.    i, k, l, m        : integer ;
  160. begin
  161.    zz :=  0.30 ;   {factor used to control decrease of condition of matrix}
  162.    for l := 1 to 5 do
  163.       begin
  164.       xx :=  zz*(3-l) ;
  165.       X  :=  exp(xx/LOG10E);  {slowly decreases conditioning}
  166.       filla ; fillb ; fillc ;
  167.       mult ; sumit ;
  168.       err[l] := sum/sqr(N) ;  {error is average absolute error per element}
  169.       if err[l] > MINERR then logerr[l] := -ln(err[l]) * LOG10E
  170.                           else logerr[l] := LOGMIN;
  171.       testerr[1] := testerr[1] + (LOGMIN - logerr[l]) ;
  172.       end ;
  173.    testerr[1] := testerr[1]/5.0 ;
  174.  
  175.    DrawAt(5,21,'#1a: 10x10 matrix       ');
  176.    f.Style := FixedDecimal;
  177.    f.Digits := 1;
  178.    z := 26;
  179.    for i := 1 to 5 do 
  180.      begin
  181.        z := z + 5;
  182.        Num2Str(f,logerr[i],s);
  183.        DrawAt(z,21,s);
  184.      end;
  185.    f.Digits := 2;
  186.    Num2Str(f,testerr[1],s);
  187.    DrawAt(59,21,s);
  188.  
  189. { infinite product for 1-x: run in reverse to test division }
  190.  
  191.    sum := 0.0 ;
  192.    for l := 1 to 5 do
  193.       begin
  194.       xx := (1 - l) / 4.0 ;
  195.       zz := exp((xx-2.0)/LOG10E); {increases number of factors for convergence}
  196.       xx := 1.0 - zz ;  { zz ≈ .01 => loss of 2 sig figs }
  197.    {The following formula for the number of factors is designed to give
  198.       sufficient accuracy, while avoiding underflow in the powers of xx.
  199.       It gives a more uniform computation from compiler to compiler.}
  200.       m := 12+l ;
  201.       quot := 1.0 ;
  202.       for k := 1 to m do
  203.          begin
  204.          quot := quot / (1.0 + xx) ;
  205.          xx := xx * xx ;
  206.          end ;
  207.       err[l] := abs(1.0 - quot/zz)*0.01 ;
  208.       { factor of 0.01 to compensate for cancellation errors }
  209.       if err[l] > MINERR then diverr[l] := -ln(err[l]) * LOG10E
  210.                          else diverr[l] := LOGMIN;
  211.       sum := sum + (LOGMIN - diverr[l]) ;
  212.       logerr[l] := diverr[l] ; {needed for later average}
  213.       end ;
  214.    sum := sum / 5.0 ;
  215.  
  216.  
  217.  DrawAt(5,22,'#1 : infinite product   ');
  218.    f.Digits := 1;
  219.    z := 26;
  220.    for i := 1 to 5 do 
  221.      begin
  222.        z := z + 5;
  223.        Num2Str(f,diverr[i],s);
  224.        DrawAt(z,22,s);
  225.      end;
  226.    f.Digits := 2;
  227.    Num2Str(f,sum,s);
  228.    DrawAt(59,22,s);
  229.  
  230.  
  231. { test continued fraction for tangent against exact values for five angles:
  232.    this is a test of division and subtraction, not of the tangent.}
  233.  
  234.    th[1] := PI/12.0 ;
  235.    th[2] := PI/6.0 ;
  236.    th[3] := PI/4.0 ;
  237.    th[4] := PI/3.0 ;
  238.    th[5] := 5.0*PI/12.0 ;
  239.    val[1] := 2.0 - ROOT3 ;
  240.    val[2] := 1.0 / ROOT3 ;
  241.    val[3] := 1.00 ;
  242.    val[4] := ROOT3 ;
  243.    val[5] := 2.0 + ROOT3 ;
  244.    sum := 0.0 ;
  245.    m := 8 ;   { this number of iterations gives sufficient accuracy }
  246.    for l := 1 to 5 do
  247.       begin
  248.       a0 := 2.0 * m + 1.0 ;
  249.       p2 := th[l] ;
  250.       p  := sqr(p2) ;
  251.       d0 := a0 - p / (a0 + 2.0) ;
  252.       for k := 1 to m do
  253.          begin
  254.          a1 := a0 - 2.0 ;
  255.          d1 := a1 - p / d0 ;
  256.          a0 := a1 ;
  257.          d0 := d1 ;
  258.          end ;
  259.       frac := p2 / d0 ;
  260.       funct[l] := frac ;
  261.       end ;
  262.  
  263.    for l := 1 to 5 do
  264.       begin
  265.       err[l] := abs(1.0 - val[l]/funct[l]) ;
  266.       if err[l] > MINERR then diverr[l] := -ln(err[l]) * LOG10E
  267.                           else diverr[l] := LOGMIN;
  268.       sum := sum + (LOGMIN - diverr[l]);
  269.       end ;
  270.    sum := sum / 5.0 ;
  271.  
  272.  
  273.    DrawAt(5,23,'#1 : continued fraction ');
  274.    f.Digits := 1;
  275.    z := 26;
  276.    for i := 1 to 5 do 
  277.      begin
  278.        z := z + 5;
  279.        Num2Str(f,diverr[i],s);
  280.        DrawAt(z,23,s);
  281.      end;
  282.    f.Digits := 2;
  283.    Num2Str(f,sum,s);
  284.    DrawAt(59,23,s);
  285.    
  286.    for i := 1 to 5 do
  287.       begin
  288.       logerr[i] := 0.5*(logerr[i] + diverr[i]) ;
  289.       testerr[2] := testerr[2] + (LOGMIN - logerr[i]);
  290.       end ;
  291.    testerr[2] := testerr[2]/5.0 ;
  292.  
  293. DrawAt(5,24,'#1b: division average   ');
  294.    f.Digits := 1;
  295.    z := 26;
  296.    for i := 1 to 5 do 
  297.      begin
  298.        z := z + 5;
  299.        Num2Str(f,logerr[i],s);
  300.        DrawAt(z,24,s);
  301.      end;
  302.    f.Digits := 2;
  303.    Num2Str(f,testerr[2],s);
  304.    DrawAt(59,24,s);
  305.  
  306. end;
  307.  
  308.  
  309. procedure trig ;   {TEST 2: first, errors in some sine identities }
  310. var
  311.    i, j, k, l        : integer ;
  312. begin
  313.    for l := 1 to 5 do logerr[l] := 0.0 ;
  314.    for j := 1 to ITERTRIG do
  315.       begin
  316.       p := j - 1 ;
  317.       th[1] := PI/12.0 + p*PI ;
  318.       th[2] := PI/6.0 + p*PI ;
  319.       th[3] := PI/4.0 + p*PI ;
  320.       th[4] := PI/3.0 + p*PI ;
  321.       th[5] := 5.0*PI/12.0 + p*PI ;
  322.       val[1] := osgn(j-1)*ROOT2*(ROOT3-1.0)*0.25 ;
  323.       val[2] := osgn(j-1)*0.5 ;
  324.       val[3] := osgn(j-1)*SQRTO2 ;
  325.       val[4] := osgn(j-1)*0.5*ROOT3 ;
  326.       val[5] := osgn(j-1)*ROOT2*(ROOT3+1.0)*0.25 ;
  327.       for l := 1 to 5 do funct[l] := sin(th[l]) ;
  328.       for l := 1 to 5 do
  329.          begin
  330.          err[l] := abs(1.0 - val[l]/funct[l]) ;
  331.          if err[l] > MINERR then logerr[l] := logerr[l] - ln(err[l]) * LOG10E
  332.                             else logerr[l] := logerr[l] + LOGMIN ;
  333.          end;
  334.       end;
  335.    for l := 1 to 5 do logerr[l] := logerr[l] / ITERTRIG ;
  336.    for l := 1 to 5 do testerr[3] := testerr[3] + (LOGMIN - logerr[l]);
  337.    testerr[3] := testerr[3]/5.0 ;
  338.  
  339. DrawAt(5,25,'#2a: sine function      ');
  340.    f.Digits := 1;
  341.    z := 26;
  342.    for i := 1 to 5 do 
  343.      begin
  344.        z := z + 5;
  345.        Num2Str(f,logerr[i],s);
  346.        DrawAt(z,25,s);
  347.      end;
  348.    f.Digits := 2;
  349.    Num2Str(f,testerr[3],s);
  350.    DrawAt(59,25,s);
  351.  
  352. { compare sin()/cos() with exact values for tan() }
  353.    for l := 1 to 5 do logerr[l] := 0.0 ;
  354.    for j := 1 to ITERTRIG do
  355.       begin
  356.       p := j - 1 ;
  357.       th[1] := PI / 12.0 + (j-1)*PI ;
  358.       th[2] := PI / 6.0 + (j-1)*PI ;
  359.       th[3] := PI / 4.0 + (j-1)*PI ;
  360.       th[4] := PI / 3.0 + (j-1)*PI ;
  361.       th[5] := 5.0 * PI / 12.0 + (j-1)*PI ;
  362.       val[1] := 2.0 - ROOT3 ;
  363.       val[2] := 1.0 / ROOT3 ;
  364.       val[3] := 1.00 ;
  365.       val[4] := ROOT3 ;
  366.       val[5] := 2.0 + ROOT3 ;
  367.       for l := 1 to 5 do funct[l] := sin(th[l])/cos(th[l]) ;
  368.       for l := 1 to 5 do
  369.          begin
  370.          err[l] := abs(1.0 - val[l]/funct[l]) ;
  371.          if err[l] > MINERR then logerr[l] := logerr[l] - ln(err[l]) * LOG10E
  372.                             else logerr[l] := logerr[l] + LOGMIN ;
  373.          end  ;
  374.       end;
  375.    for l := 1 to 5 do logerr[l] := logerr[l] / ITERTRIG ;
  376.    for l := 1 to 5 do
  377.       testerr[4] := testerr[4] + (LOGMIN - logerr[l]) ;
  378.    testerr[4] := testerr[4]/5.0 ;
  379.  
  380. DrawAt(5,26,'#2b: tangent|sin()/cos()');
  381.    f.Digits := 1;
  382.    z := 26;
  383.    for i := 1 to 5 do 
  384.      begin
  385.        z := z + 5;
  386.        Num2Str(f,logerr[i],s);
  387.        DrawAt(z,26,s);
  388.      end;
  389.    f.Digits := 2;
  390.    Num2Str(f,testerr[4],s);
  391.    DrawAt(59,26,s);
  392.  
  393.  
  394. { compare arctan() with tan() for consistency }
  395.    for l := 1 to 5 do logerr[l] := 0.0 ;
  396.    for j := 1 to ITER do
  397.       begin
  398.       for l := 1 to 5 do th[l] := (5*j+l-5)*PIO2/(5*ITER+1) ;
  399.       for l := 1 to 5 do val[l] := sin(th[l])/cos(th[l]) ;
  400.       for l := 1 to 5 do funct[l] := arctan(val[l]) ;
  401.       for l := 1 to 5 do
  402.          begin
  403.          err[l] := abs(1.0 - th[l]/funct[l]) ;
  404.          if err[l] > MINERR then logerr[l] := logerr[l] - ln(err[l]) * LOG10E
  405.                             else logerr[l] := logerr[l] + LOGMIN ;
  406.          end;
  407.       end;
  408.    for l := 1 to 5 do logerr[l] := logerr[l] / ITER ;
  409.    for l := 1 to 5 do
  410.       testerr[5] := testerr[5] + (LOGMIN - logerr[l]);
  411.    testerr[5] := testerr[5]/5.0 ;
  412.  
  413. DrawAt(5,27,'#2c: arctan function    ');
  414.    f.Digits := 1;
  415.    z := 26;
  416.    for i := 1 to 5 do 
  417.      begin
  418.        z := z + 5;
  419.        Num2Str(f,logerr[i],s);
  420.        DrawAt(z,27,s);
  421.      end;
  422.    f.Digits := 2;
  423.    Num2Str(f,testerr[5],s);
  424.    DrawAt(59,27,s);
  425.  
  426. end;
  427.  
  428.  
  429. procedure transc ;   {TEST 3: ln() and exp() for consistency}
  430. var
  431.    i, j, l        : integer ;
  432. begin
  433.    for l := 1 to 5 do logerr[l] := 0.0 ;
  434.    for j := 1 to ITER do
  435.       begin
  436.       for l := 1 to 5 do th[l] := (5*j+l-5)*STEP ;
  437.       for l := 1 to 5 do val[l] := exp(th[l]) ;
  438.       for l := 1 to 5 do funct[l] := ln(val[l]) ;
  439.  
  440.       for l := 1 to 5 do
  441.          begin
  442.          err[l] := abs(1.0 - th[l]/funct[l]) ;
  443.          if err[l] > MINERR then logerr[l] := logerr[l] - ln(err[l]) * LOG10E
  444.                             else logerr[l] := logerr[l] + LOGMIN;
  445.          end;
  446.       end;
  447.    for l := 1 to 5 do logerr[l] := logerr[l] / ITER ;
  448.    for l := 1 to 5 do
  449.       testerr[6] := testerr[6] + (LOGMIN - logerr[l]);
  450.    testerr[6] := testerr[6]/5.0 ;
  451.  
  452. DrawAt(5,28,'#3a: ln() & exp()       ');
  453.    f.Digits := 1;
  454.    z := 26;
  455.    for i := 1 to 5 do 
  456.      begin
  457.        z := z + 5;
  458.        Num2Str(f,logerr[i],s);
  459.        DrawAt(z,28,s);
  460.      end;
  461.    f.Digits := 2;
  462.    Num2Str(f,testerr[6],s);
  463.    DrawAt(59,28,s);
  464.  
  465. end;
  466.  
  467.  
  468. procedure roots ;  { sqrt() and sqr() identities }
  469. var
  470.    i, j, l        : integer ;
  471. begin
  472.    for l := 1 to 5 do logerr[l] := 0.0 ;
  473.    for j := 1 to ITER do
  474.       begin
  475.       for l := 1 to 5 do th[l] := (5*j+l-5)*STEP ;
  476.       for l := 1 to 5 do val[l] := sqrt(th[l]) ;
  477.       for l := 1 to 5 do funct[l] := sqr(val[l]) ;
  478.  
  479.    for l := 1 to 5 do
  480.       begin
  481.       err[l] := abs(1.0 - th[l]/funct[l]) ;
  482.       if err[l] > MINERR then logerr[l] := logerr[l] - ln(err[l]) * LOG10E
  483.                          else logerr[l] := logerr[l] + LOGMIN;
  484.       end;
  485.    end;
  486.    for l := 1 to 5 do logerr[l] := logerr[l] / ITER ;
  487.    for l := 1 to 5 do
  488.       testerr[7] := testerr[7] + (LOGMIN - logerr[l]) ;
  489.    testerr[7] := testerr[7]/5.0;
  490.  
  491. DrawAt(5,29,'#3b: sqrt() & sqr()      ');
  492.    f.Digits := 1;
  493.    z := 26;
  494.    for i := 1 to 5 do 
  495.      begin
  496.        z := z + 5;
  497.        Num2Str(f,logerr[i],s);
  498.        DrawAt(z,29,s);
  499.      end;
  500.    f.Digits := 2;
  501.    Num2Str(f,testerr[7],s);
  502.    DrawAt(59,29,s);
  503.  
  504. end;
  505.  
  506.  
  507. procedure Initialize;
  508. { Perform one-time setup work. }
  509. var
  510.   r: Rect;
  511.   MyNewWindow: WindowPtr;
  512. begin
  513.   InitGraf(@thePort);           { initialize QuickDraw }
  514.   InitFonts;                    {     "      Font Manager }
  515.   InitWindows;                  {     "      Window Manager }             
  516.   InitCursor;                   { make cursor an arrow }
  517.   done := False;
  518.   SetRect(r,30,60,600,400);
  519.   MyNewWindow:=NewWindow(nil,r,'Accuracy Tester',True,documentProc,Pointer(-1),True,0);
  520.   MyNewWindow^.txFace := [bold]; MyNewWindow^.txFont := SystemFont;
  521.   windowPeek(MyNewWindow)^.refCon := Longint(NewString('Contents of window 1'));
  522.   
  523.  end;
  524.  
  525.  procedure Header ;
  526. begin
  527.    NewFont ( monaco, 9 );
  528.    TextFace( [] );
  529.    DrawAt(5,6,'ACCURACY: extended tester: ' + COMPIL + '; ' + MACHIN) ;
  530.    DrawAt(5,7,'  V 1.7 (c) 1987, PC Tech Journal and Ziff Communications Co.');
  531.    DrawAt(5,8,'                  written by Jim Roberts.');
  532.    DrawAt(5,9,'    Modified for Macintosh & Turbo Pascal by Morry Hodges, CIS 74766,3426');
  533.    DrawAt(5,10,'Test 1 checks multiplication and addition, then division and subtraction.');
  534.    DrawAt(5,11,'Test 2 measures the accuracy of the Turbo trig functions (there is no tan()).');
  535.    DrawAt(5,12,'Test 3 finds the truncation error in some exponential and sqrt identities.');
  536.    DrawAt(5,13,'ACCURACY is the rounded negative log of error.  Program may exit abnormally.');
  537.    DrawAt(5,14,'NOTE: an increase of 1 in the rating means a factor of TEN less accurate.');
  538.    DrawAt(5,16,'Interpretation  <0.0 - 0.5 => Excellent     1.0 - 1.5 => Fair');
  539.    DrawAt(5,17,'  of RATING:     0.5 - 1.0 => Good          1.5 <     => Poor');
  540.    DrawAt(5,18,'') ;
  541.    DrawAt(5,19,'      TESTS                       ACCURACY           RATING ');
  542. end;
  543.  
  544. begin  { program ACCURACY begins here }
  545.   Repeat
  546.    Initialize;
  547.    Header ;
  548.    for i := 1 to  10 do testerr[i] := 0.0 ;
  549.    arith ;
  550.    trig ;
  551.    transc ;
  552.    roots ;
  553.    ntest := 7 ;
  554.    toterr := 0.0;
  555.    for i := 1 to ntest do toterr := toterr + testerr[i] ;
  556.    toterr := toterr/ntest ;
  557.    DrawAt(5,31,'Overall rating:');
  558.    Num2Str(f,toterr,s);
  559.    DrawAt(23,31,s);
  560.    DrawAt(5,33,'Click mouse to finish');
  561.    Pause;
  562.       until Button;
  563. END.
  564.  
  565.  
  566.